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ABSTRACT 

We present a new collocation method for the numerical solution of partial differential 
equations. This method uses the Chebyshev collocation points, but because of the way the 
boundary conditions are implemented, has all the advantages of the Legendre methods. In 
particular L 2 estimates can be obtained easily for hyperbolic and parabolic problems. 
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1 Introduction 


Polynomial pseudospectral (or collocation) methods have been extensively used in the nu- 
merical solutions of partial differential equations. The underlying idea in those methods is 
to approximate the unknown function by an interpolation polynomial at some pre-described 
(collocation) points. The polynomial is then required to satisfy the PDE at the collocation 
points. This procedure yields a system of ordinary differential equation to be solved. 

Historically, (see [10]) the first such points to be used were the Chebyshev collocation 
points 

Xj = cos(^) 0<j<N. 

Those points were chosen because they allowed the use of Fast-Fourier-Transforms in the 
computations. It was only later (see [7]) that those points were identified with the nodes 
of the Gauss- Lobatto- Chebyshev (G-L-G) quadrature formula. This observation is the key 
in the stability analysis of the pseudospectral Chebyshev methods. The G-L-C quadrature 
formula led to the weighted L 2 norm 



However, it has been noted in [8] that this is not a natural norm for hyperbolic equations. 
In fact the differential equation is not well posed in this norm. Also it complicated the 
stability analysis even for parabolic equations. The theory (and therefore the confidence in 
applying those methods) is not complete. 

Once the connection between the collocation points and the Gauss Lobatto points is 
established, it is natural to use the nodes of the Gauss- Lobatto - Legendre (G-L-L) quadrature 
formula. We refer the reader to [2] for review of those methods. Recently [1] an 0(N log N) 
method was proposed for the Legendre points. The main problem with those points are that 
they are not given explicitly, and their evaluation for large N is not robust due to roundoff 
errors. 

In this paper we present a method (and name it The Chebyshev- Legendre Method) that 
has the advantages of both the Chebyshev and Legendre methods. The method utilizes the 
Chebyshev collocation points allowing the use of fast Fourier algorithms and avoiding the 
roundoff error associated with computing the Legendre grid points. The boundary conditions 
are imposed via a new penalty technique in such a way that the method is stable in the usual 
Li norm (rather than the weighted L 2 norm). Hence the Chebyshev- Legendre method enjoys 
the advantages of the Chebyshev method as well as those of the Legendre method. 
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The implementation of the boundary conditions is done by a penalty method. A penalty 
term is added to the PDE at all grid points in such a way that, in the limit of number of 
grid points tend to infinity, the boundary conditions are satisfied. This procedure seems to 
be better than the direct imposition of the boundary conditions, and in our case has the 
extra advantage of yielding the Legendre method at the Chebyshev points. 

A similar idea had been tried by Reyna [11], The difference between his approach and 
ours is in the imposition of the boundary conditions. Instead of transforming from the 
Chebyshev basis to the Legendre one as in [11], we impose the boundary conditions via 
penalty method and through that automatically switch to the Legendre basis without using 
it in the differentiation procedure. 

The paper is organized as follows: 

In Section 2 we quote the essential formulas for the use of Chebyshev and Legendre 
methods. 

In Section 3 we present the Chebyshev-Legendre method for hyperbolic equations. In 
subsection 3.1 we describe the method and prove (in Theorem (3.1.1)) an energy estimate 
to show stability. In Theorem (3.1.2) we bring another version of this energy estimate. 

In subsection 3.2 we consider the relationship of the new method to the Legendre penalty 
method and show that the differentiation matrices of the two methods are related via a 
similarity transformation. This fact is proved in Theorem 3.2.2 

In Section 4 we discuss and prove the stability of Chebyshev-Legendre method for the 
heat equation with Robin boundary conditions. 

Section 5 concludes the paper with some numerical experimentations with the new 
method. 

In future work we will report on the convergence results of the new method for nonlinear 
hyperbolic equations. 


2 Preliminaries 


This Section is devoted to the definitions of the pseudospectral methods to be used later. We 
will discuss Chebyshev and Legendre methods which are based on the Chebyshev polynomials 

T n (x) = cos(A^cos _1 x) (2-1) 


and the Legendre polynomials 


P N (x) 


1 d N 
2 N N\ dx N 


(x'-l)" 


( 2 . 2 ) 


respectively. Associated with these two polynomials are several Gauss-type quadrature for- 
mulas. We will consider, in this paper, the Gauss Lobatto type formulas. 
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We start by defining the Chebyshev Collocation points Xj by 


x i = cos (-jf) 0<j<N 


(2.3) 


These points are the zeroes of the polynomial (1 — x 2 )T' N (x) and associated with it, we have 

The Gauss Lobatto Chebyshev Quadrature Formula: 

Let f(x) be a polynomial of degree 2 N — 1, then 


3 = 0 1 


(2.4) 


where the weight Cj are given by 


c i = 


Co — C N 


7T 

7v 

IT 

2 77 


l <j < N 


(2.5) 


Similarly, the Legendre Collocation points yj are defined as the roots of the polynomial 
(1 — x 2 )P' n (x). For these points we have 

The Gauss Lobatto Legendre Quadrature Formula: 

Let f(x) be a polynomial of degree 2 N — 1, then 


E/( yj)^- f 


(2.6) 


J=0 

where the Gauss Lobatto weights ojj are given by 

2 


yv+ i 

2 




(2.7) 


u 0 = 0J N = 


N(N+ 1) 


Unlike the Chebyshev points that are known explicitly, there is no explicit formula for the 
Legendre points y 2 , they have to be computed numerically. It is interesting though that there 
is a simple formula, easily and robustly computed, for the values of the Legendre polynomials 
and their derivative at the Chebyshev points. In fact we have the following explicit formula 
for Ph(xj), (taken from [4], page 180) 


N - 1 1 3 ' 




( 2 . 8 ) 
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Pseudospectral (or Collocation) methods are based on interpolations at the points Xj or 
yj. Consider the polynomials 

Ql(x) = (l-x 2 )Pf,,(x) 

Qc(x) = (1 - x 2 )T r N (x ) , 


and define the Legendre- Lagrange polynomials by 

Ql(x) 


hj(x) 


(x - y } )Q' L {yj) ’ 
and the Chebyshev-Lagrange polynomials by 

/. x Qcjx) 

9AI (x-x,)Q' c (x,)- 

Tlien, the Legendre interpolation operator Ii is defined by 

(hf)(x) = £/(»)*,(*) 

;=o 

whereas the Chebyshev interpolation operator Ic is defined by 

TV 


By definition, we have 


j=o 


(a/)(») = fivj) 

Vcf)(Xj) = /( Xj). 


(2.9) 


( 2 . 10 ) 


( 2 . 11 ) 


( 2 . 12 ) 


From the definition of the interpolation operators Ji and Ic we get the spectral differen- 
tiation matrices T>i and T>c as follows: 

The Pseudospectral Legendre Differentiation Matrix is defined by 

(®l)jjb = KM ■ (2.13) 


The Pseudospectral Chebyshev Differentiation Matrix T>c , is defined by 


tPc)i* - g'k(x,). 


(2.14) 


(see [3] for explicit expressions for the matrices) 
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3 Hyperbolic Equations 

3.1 Scalar Hyperbolic Equation 

In this Section we consider the scalar initial-boundary value hyperbolic equation 

Ut = U x - 1 < x < 1 t> 0 (3.1) 

with the initial condition 

(/(x,0) = /(x) (3.2) 

and the boundary condition 

U(l,t)=g(t) . (3.3) 

The Chebyshev-Collocation (Pseudospectral) method involves seeking an Nth degree 
x-polynomial u N (x,t ) that satisfies 

= du "( x 4 l at x = Xj I < j < N (3.4) 

dt dx 

with the boundary condition 

u N (\,t) = u N (x 0 ,t) = g{t) (3.5) 

where x, are determined in (2.3) (xo = 1). 

Note that the equation is satisfied at all the grid points except at the boundary point 
x — 1 where the boundary condition is satisfied. 

In general, the term 4- u N (x,t ) is evaluated at all the grid points, with the use of either 
FFT or matrix-vector multiplication using the matrix T>c- Equation (3.4) is then advanced 
at all the grid points. The value of the solution at the boundary is then updated using (3.5). 

In [5] a penalty type method was introduced. In that approach we still use equation 
(3.4) for the inner points Xj, 1 < j < N, however instead of using (3.5) for the boundary, 
the following equation is satisfied 

■^— L = — ^ U =1 -r(u N (l,t) - g{t)) (3.0) 

where r is determined from stability considerations. In particular it had been found that 
stability follows if 


Equations (3.4) and (3.6) can be combined into a single equation by noting that the 
collocation points Xj defined in (2.3) are the zeroes of the polynomial (1 — x 2 )T' N (x). Thus 
the penalty method [5] can be written now as 


du N (xj,t) du N (x,t ) 


dt 


Ox 


(1 + Xj ;)T' N (xj) 

2^(1) 


(u N (\,t.) - g(t)) 


(3.7) 


for j = 0, . . . , N. 

The main difference between the penalty method (3.7) and the usual Chebyshev method 
given in (3.4) and (3.5) is that the numerical solution u w (x,t) does not satisfy the boundary 
condition exactly, but only in the limit as N ^ oo. The boundary condition is now part of 
the equation. 

Another penalty method, based on the Legendre points y 3 is presented in [6]. Similar to 
(2.6) we write this method as 


du N (y } ,t) du N {x,t) _(1 + y 3 )P' N (yj), „ ^ 

— “ “aT -1 ”" - r 2 w) M,) ' ,( 0> 

for j ~ 0, . . . , N. 

The parameter r is determined by the stability requirement. Thus the differential equa- 
tion is satisfied at the points yj,j= 1 , . . . , N. At the boundary xq = \ one uses a combination 
of the boundary condition and the differential equation. 

An obvious disadvantage of the method in (3.8) is that it utilizes the Legendre points. 
However comparing (3.7) and (3.8) shows us how to utilize the Legendre penalty method 
(3.8) at the Chebyshev points. 


The Chebyshev-Legendre (C-L) Method 

Let Pn(x ) be the Legendre polynomial of degree N. In the C-L method we seek a 
polynomial of degree N in x that satisfies 




dt 


dx 


r—'r. T 


(1 d~ • c j)T)y( j: j) 

2 ^( 1 ) 


(u N (\,t) - g(t)) 


(3.9) 


for .7 = 0,..., N. 

Note that the penalty term - — 2 p' (i) ' ’ different from zero for all the Chebyshev grid 
points Xj. Note also that applying (3.9) entails the use of the differentiation matrix T>c at the 
Chebyshev points. In fact, given u N (xj,t ) one finds the derivative based on the Chebyshev 
points, and then add the penalty term with different weights at every grid point. The term 
Pn(xj) is evaluated using the explicit formula (2.8). This is done once and for all for any 
grid size N. 
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The surprising fact is that the C-L method, though computed at the Chebyshev points, 
is stable in the usual L 2 norm, rather than the weighted L 2 norm. In fact one can state 


Theorem 3.1.1 : The h 2 Stability of the C-L Method. 

let u N {x,t ) be the solution of (3.9). Let u>j be the weights of the Gauss Lobatto Legendre 
quadrature formula and yj the nodes of the same quadrature formula. Let g(t) = 0 in (3.3) 
and (3.9), then for 


r > — = N(N + 1) 

2oJo 

the C-L method is stable in the L 2 norm. More specifically, 


N 


i= o 


N 


S *4(^,0)^- - 

j = 0 


[ [u^(l,t)(2o» 0 r - I) + u*(-l, *)]<#• 
Jo 


(310) 


Proof: 

It follows from (3.9) that 


dr^QM) = dv^jx C) _ (1 + x)P^{x) . o(f)). (3.11) 

8t dx 2P^(1) V V ' 

This is because both sides of (3.9) are polynomials of degree N that agree at N + 1 points, 
namely at the Chebyshev collocation points 0 < j < N. 

We now read (3.11) at the Legendre points y } to get 


;=0 j =0 UX 

Since the Gauss-Lobatto-Legendre quadrature formula is exact for polynomials of degree 
2N — 1 it follows that 


2 E «„(», t) ^uA Uj = 

3 = 0 

= «*(!,<) 


and thus the stability estimate (3.10) follows. The Theorem is thus proven. 


Note that unlike the Legendre-Penalty method (3.8), in which one needs to use the 
Legendre points yj in the computations, these points do not appear in the computations in 
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the C-L method. They are just introduced for the sake of the proof. The actual computations 
are done using the Chebyshev grid points xj. 

An energy estimate based on the Chebyshev points xj can be derived by using (2. 1 1) and 
(2.12) as follows: 


Theorem 3.1.2 : 
Let 


N N 


j= o 1=0 


with 


N 


= J2 9j(yk)gi(yk)^k ■ 


k = 0 


where the Chebyshev- Lagrange polynomials gj(x) are defined in (2.10). 
Then u N (-,t) satisfies the energy estimate 


(3.12) 


IM-,<)ll J = IM-.o)|| 2 - |V„(i,f)(w-i) + <(-!,()}* ( 3 . 13 ) 


Proof : 

Equation (3.13) is really a restatement of (3.10). Since u N (x,t ) is a polynomial of degree 
N in x, it can be represented exactly by 

N 

u N (x, t.) = l c u N = Y u A x h t)gi( x ) ■ 

1 = 0 

Thus 

N 

u N (yj,t) = E%K f )S'fc)’ (3-14) 

1=0 

The estimate (3.13) follows from (3.10) upon substituting (3.14) for the values of u N (t/j,<). 
The theorem is proven. 

The C-L method can be viewed from many different points of view. Theorem (3.1.1) 
shows that this method, in the constant coefficient case, is equivalent to the Legendre - 
penalty method introduced in [5]. 

Thus the C-L method is the realization of the Legendre method at the Cheby- 
shev points. 

We close this section by pointing out that the estimate (3.13) enables one to pass to 
systems easily. Actually it had been done in [5]. One has to follow the same steps. It shows 
the C-L method is stable for constant coefficients systems of hyperbolic equations. 
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3.2 The Differentiation Matrices 


Perhaps more insight can be gained if one compares the differentiation matrix induced by 
the C-L method Vcl with the differentiation matrices induced by the Chebyshev penalty 
method (3.7) T>c , and the Legendre penalty method (3.8) T>l,. 

We start by noting that the differentiation matrices T>i and T>c defined in (2.13) and 
(2.14) do not take into account any boundary conditions. The differentiation matrices in- 
duced by (3.7), (3.8) and (3.9) are variations of the basic matrices T>i and Vc, differing only 
in the method of imposing boundary condition. 

Not surprisingly, 'Dl and T>c are similar, after all both differentiate exactly polynomials 
of degree N. Since for these polynomials, the operators / l and Ic are the same, the matrices 
Vi and T>c represent the same operation in a different basis. This implies a similarity 
relationship. More specifically, this relationship can be written explicitly: 

Theorem 3.2.1 : 

Let S be the matrix whose elements Sj t k are given by 

s**-M**) (s-w) 

where hj( x) are the Legendre- Lagrange polynomials defined in (2.9). Again Xj are the 
Chebyshev points. 

Let T be the matrix whose elements Tj t k are given by 

Tj,k — ( 3 - 16 ) 

where g 3 (x) the Chebyshev- Lagrange polynomials are defined in (2.10). i/a., are the Legendre 
collocation points. 

Then 


S = T ~ 1 (3.17) 

and 

V c = SV l T (3.18) 


Proof : 

Since Qj{x) is a polynomial of degree N it is given by 

9ii x ) = J29j(yi) h i( x )- ( 3 - 19 ) 

(=0 
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Substituting x*. and making use of the fact that g 3 {xk) = ftj,k , we get 

N 

f>j,k = 

proving that 

I = ST . 

Differentiating (3.19) we get 

afa) = I Zsj(yt) h i(x) 

1=0 

However /tj(x) is itself a polynomial of degree N and therefore it can be expressed as 

h\(x) = h 'i{ym)h m (x) , 

m=0 

which leads to 

g'jM = 53 S 9i(yi) h 'i(y™)K{xk) 

1—0 TH — 0 

The theorem is thus proved. 

We will show now that the differentiation matrix induced by the C-L method T>cl I s 
similar to the differentiation matrix T>cp induced by the Legendre penalty method (3.8). 
This will demonstrate the fact that the C-L method is the realization of the Legendre method 
on the Chebyshev grid. 

Theorem 2.2.2: 

Let T>cl be the differentiation matrix induced by the Chebyshev Legendre method (3.9) 
and T>lp the differentiation method induced by the Legendre penalty method (3.8) then 

V CL = SV LP J (3.20) 

where the matrices <5, 7” defined in (3.15), (3.16) are the transformation matrices between 
the Chebyshev points and the Legendre points. 

Proof : 

Note that the differentiation matrix Vpp is essentially the matrix Vp introduced in (2.13), 
modified to take into account the boundary conditions, imposed via penalty in (3.8). Thus 

(PLp)j,k = ~ xSo^j.o ■ (-^ 21 ) 
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In the same manner we can write explicitly 


m \ (r > ^ (1 +x k )F ) ' N (x k ) 

{VcL)j,k = c)j,k ~ r 2Py(l) d ° ' 


(3.22) 


Equation (3.22) is a direct consequence of (3.9). Note that the full first column of T>c is 
modified, and not only the first element as in (3.21). 

We proceed by writing explicitly the elements of the matrix ST^i,f‘T . In fact 


N N 


(ST>Lf>T)j y k — {ST>lT) h k ~ T 9j{yi)^0,lfio,mh, n { x k) ■ 

1=0 m= 0 


Thus using (3.18) we get 


{ST> L pT)j tk = {V c )j,k - rgjiyo^oixk) 


(3.23) 


From (2.9) 


h 0 (Xk) = 


(1 + Xk)P]y(Xk) 


2^(0 

and since y Q = x 0 = 1 g j(y 0 ) = Sj,o so the right hand side of (3.23) is exactly the same as 
this of (3.22). 

Thus (3.20) is established. The proof is completed. 


4 Parabolic Equations 


In this Section we present the Chebyshev- Legendre method for the parabolic equation 

du d 2 u 


3t dx 2 

with Robin boundary condition 


1 < x < 1 ,t > 0 


(4.1) 


+ = g + (t) 

7u(-l,<) + 6u s (-l,i) =<7 _ (f). 


(4.2) 


We will assume that o',/?, 7 are non-negative and 6 is non-positive. This assures the time 
decay (or non-growth) of u(x,t). 

We note that by now there is a very limited stability theory for the Chebyshev method. 
In fact stability had been proved first for the Dirichlet case f) = 0, 5 — 0, (see [7]) and then 
for Neumann case 0 = 7 = 0. [9]. Here we present the C-L method and prove stability for 
the approximation to (4.1), (4.2) for the general Robin case. 
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Denote by Vn the finite dimensional space of polynomial of degree at most N. We define 
the operator A 


A : Vn — > Vn 


by 

Av(x,t) = -^A + *(*.<) («) 

where 


R(x,t) = r 0 g+(x)[B + (i) - 9 + (()] + r„Q-(x)\B-(t) - g~(t )] (4.4) 


with 

<?+(*) = ( ‘ ’ B+W = HM) + /MM) 

<3~M = (1 wl(i) X) ' B ' {t)= 7»(-l ,0 + M-M). 

The numbers r 0 , r N will be determine later to assure stability. We define also the following 
scalar product 


N 

(v, w)n = Y v (Vj) w (yi) u j ( 4 - 5 ) 

j=0 

where yj,u>j are the Legendre points and weights respectively. 


The Chebyshev-Legendre Method for Parabolic Equations 

We seek the polynomial of degree N in x, v(x,t) that satisfies 
dv(xj,t ) d 2 v(x,t). 


dt dx 2 11-1 

where x } are Chebyshev collocation points. 


R(xj,t) 0 <j<N 


(4.6) 


Note that again, the work is done on the Chebyshev points xj, the penalty values 
<5+(xj), are computed by (2.8) and are nonzero for any Xj. 

To prove the stability of (4.6), we set g + (t.) = g~(t) = 0. In the following lemma, we will 
find conditions on r 0 and t n such that the operator A to be semi bounded. 

Lemma 4.1 : 
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Let v £ Vn and 


T a,b — ^-^[(1 + 2 k ) + 2 \/k + K 2 ] 
T a,b = ^-^[(1 + 2 k) — 2 yfn + k' 2 ] 

with k = u>o(i/b. 

Let the operator A be defined in (4.3). Then 


N - 1 


{Av,v) n > v x(yj)uj 
j=i 


provided 


T *,P < T 0 < T+p 

r ;,\6\ < T » <K\S\‘ 


(4.7) 


(4.8) 

(4.9) 


Proof : 

Since the Gauss Lobatto Quadrature formula (2.6) is exact for polynomials of degree 
2 N — 1 and since v(x,t) is a polynomial of degree N, we have 


f 1 

= / v(x)v xx (x)dx 

j=o J ~ 1 

= v(\)v x (l) — v(—l)v x (—l) — J v x (x)v x (x)dx 

using the standard integration by parts technique. 

Using again the Gauss Lobatto formula, one would get 

N N 

-'52 v (yj) v z*(yj) u j = - v ( 1 K( 1 ) + v (-i) w *(-i) 

j=0 J-O 

TV— 1 

= J2 v ' 2 MPi + 

J=i 

^(l)wo + U*(-l )u N - u(l)n x (l) + u(-l)Ua:(-l) 

Thus making use of (4.10) we can write 


(4.10) 


{Av,v) N = F(l,a, 0) + F(-\ N) + ]T) v 2 x {y } )uj 


j = i 


(4.11) 
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where 

F(x,a,b,k) = x(r k bu> k - l)u(x)u x (x) + T k au k v 2 (x) + w k v 2 (x) . (4-12) 

In order for A to be positive we need to choose r 0 and r N such that F(l , a, ft, 0) and 
F(-I,7, |6|, N) are non-negative. For F(l,a,(3,0) to be positive, we need 

( TofiuJo ~ l) 2 < 4o r 0 a.’o 
or 

t 2 0W q — 2toWo(/^ + 2cvu>o) + 1 < 0 

Thus T 0 has to lie between the roots of the parabola described in the left hand side, 
namely t~ 0 and 

The same kind of consideration holds for r N . Thus F(l,cv,/^,0) and F(— 1 , 7 , |6|, N) are 
non-negative for the range of r 0 and r N given in (4.8) and (4.9), respectively. (4.7) follows 
from (4.11). 

Remarks 

1. The Dirichlet boundary condition for x — 1 is obtained from (4.2) by setting a = 
l,/4 = 0. In this case 

r ito = 00 
^ = 

which yields the condition for the penalty amplitude 

r 0 > ^N 4 (N+l) 4 

2. The Neumann boundary condition for x = 1 corresponds to the case a - 0,/3 = 1. In 
this case r 0 + — Tq X yielding the condition 

1 N(N+\) 

To = —= 7 

LOq l 

We are now ready to state the stability theorem for the C-L method when applied to 
parabolic equations with Robin boundary conditions : 

Theorem 4.1 : 
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Let r 0 and t n satisfy (4.8) and (4.9) respectively. Let v(x,t.) 6 Vn be the C-L approxi- 
mation to u(x,t), obtained by (4.6). Assuming that <7 + (t) = = 0, v(x,t) satisfies the 

energy estimate 

rT N-1 

( v(x,t),v(x,t)) N < (v(x,0),v(x,0)) N -2 / v x{yj^') dt (4*13) 

Jo j=i 

where the scalar product (/, g)pj is defined in (3.5). 


Proof : 

Since (4.6) holds for j = 0, ..., N and since v,v xx and R are polynomials of degree at most 
N, we conclude that both sides of (4.6) are equal not only at the grid points but also for 
every x. 


dv(x,t.) d 2 v(x,t 


dt 


dx 2 


— R(x, t ) 


- 1 < x < 1 


where R(x,t.) is defined in (4.4). 

Noting the definition of A in (4.3), we get 

dv 


3t 


Av . 


Thus 


Llsing Lemma 4.1 yields 


( v , = (v,Av)n 


Id N ~' 

2 di^ V,V ^ N ~ ~ ? 


and integration yields the stability result (4.13). 


We stress again that the Legendre collocation points y 3 are ’’ghost points”, which are 
never used in the computations but only in the proof of the stability. Actually we could 
restate the proof in terms of the Chebyshev collocation points Xj as in Theorem 3.1.2. 

5 Numerical Results 

Case 1: Linear scalar PDE 
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In this Section, we will consider some numerical examples that verify our claims stated 
in previous Sections. Consider the scalar linear initial-boundary value hyperbolic PDE 

U t = U x -l<x<l,t>0 (5.1) 

with initial condition 

f/(x,0) = sin(27rfcx) 

and boundary condition at x = 1 

= g(t) = sin(27rfc(l + t )) 

We seek an N degree. x-polynomial v(x,t.) that satisfies 

^u(xj,t) = Dv(xj,t) - TQ(xj)(v(\,t) - g(t)) (5.2) 

at Chebyshev collocation point Xj = cos(tt j/N), j = 0, . . . , N and D is the differentiation 
operator (matrix). 

For different construction of the A^tli degree polynomial Q(x), one could have different 
type of boundary treatments. For examples, 

1. if Q(x) = {}+ 2 p, F (t ) l) , r> = D c and x = x 3 are the Gauss- Lobatto-Chebyshev points, 
then we have the Chebyshev- Legendre method (C-L). 

2. if Q{x) = D = Dc and x - x i are tlie Gauss- Lobatto-Chebyshev points, 

then we have the Chebyshev penalty method (C-P). 

3. if Q(x) = , D = D l and x = y 3 are the Gauss-Lobatto-Legendre points, then 

we have the Legendre penalty method (L-P). 

Let denote = v(xj,t n ) and A t be the time step increment, then for j = 0,..., N, 
we would advance the system of ODE (5.2) in time by the third order Heun Runge Kutta 
scheme that has the following form : 

For j = 0,1,..., A r , and = g(0) 


, 0 ) _ „(*) 


v 


V 


( 2 ) 
J 

(n+1) 


v r + 


= v\ n) + - rQ(xj)(v'j l> - g(t n ) - -V(/„))) 


At 


(Dv\ n) - TQ(xj)(v] n> - g(t n ))) 


3 
2A t 


,0) 


„0) 


At 


3 

1 (h) , 3 (t) 

+ 4 V ) + 

3A t 


(5.3) 


p) _ rO( r. \(J 2 ) _ „(t \ - \ - ^Lg"(t n ))) 


{Dv) - tQ{x 3 )(v) - g{t „) - — g (<„) - 


9 
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where r/(/ n ) and g r/ (t n ) are the derivatives of the time-dependent boundary conditions in 
time at t = t n . 

It has been observed before that if one imposes boundary condition at each intermediate 
stages of the Runge Kutta scheme, a larger time-step (CFL number) can be used. Otherwise, 

CFL number has to be reduced by as much as four time for stability. In this study, we define 

At = CFL/N 2 . 

The traditional way of the imposing exact boundary condition at x = 1 can be described 
as following : 

For j = 0, . . . , N y and Vj ° ^ = f/(x, 0) 

” 1 " = + y''"!" 1 

t'o I) = 

vf> = „<”> + 

»S 2) = g(K+'^) (5.4) 

(u+1) 1 (n) . 3 (1) 3 A t (2) 

v ) = 4 V ) + 4 } + ~T DV 3 
u o’ l+1) = + A t) 

However, as shown in Table I that this procedure would lead to reduction of accuracy in 
time as N increases. 

Table I 

L -2 Error and order of accuracy for (5.4) with k = 1 


N 

Error 

Rate 

Error 

Rate 

Error 

Rate 

16 

0.82E-03 


0.10E-03 


0.29E-05 


32 

0.15E-04 

2.89 

0.18E-05 

2.91 

0.28E-07 

3.35 

64 

0.42E-06 

2.57 

0.49E-07 

2.61 

0.72E-09 

2.64 

128 

0.17E-07 

2.31 

0.19E-08 

2.33 

0.28E-10 

2.34 

CFL 

8 


4 

1 


Hence, the above procedure is modified as following (Detailed discussion and analysis 
will appear in a future paper) : 

For j = 0,1,..., A^, and = U(x, 0) 



(») 

j 



(n) 

i 
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(5.5) 


.i (1 > 

u o 


g(Q + y »'(<») 


( 2 ) _ » 


V) — V 


+ 


2 At 


3 Dv > 


0) 


J2) _ 
u 0 — 

,(«+!) 

i 

,(«+!) _ 


, , 2A t , 2At 2 , 
(]{tn) H — — 9 (Li) 3 „ <7 (Li) 


1 


,00 


3 

3 , 0 ) 


4 v r + i v j + 
<7(<«+i) 


9 

3A ; 
-r D,, 7 


( 2 ) 


We shall denote this procedure as (XBC). Table I indicated that the order of time accuracy 
for this procedure is third order for all N. 


Table II 

L -2 Error and order of accuracy for (5.5) with k — 1 


N 

Error Rate 

Error Rate 

Error Rate 

16 

32 

64 

128 

0.77E-03 
0.12E-04 3.00 

0.19E-06 3.00 

0.30E-08 3.00 

0.98E-04 
0.15E-05 3.00 

0.24E-07 3.00 

0.37E-09 3.00 

0.28E-05 
0.24E-07 3.44 

0.37E-09 3.00 

0.58E- 1 1 2.99 

CFL 

8 

4 

1 


Next, using the C-L method, one get the Li error and the order of accuracy as listed in 
Table II. 


Table III 

L -2 Error and order of accuracy for C-L method with k = l,r = 4u> 0 


N 

Error Rate 

Error Rate 

Error Rate 

16 

32 

64 

128 

0.47E-03 
0.74E-05 2.99 

0.12E-06 3.00 

0.18E-08 3.00 

0.60E-04 
0.93E-06 3.01 

0.15E-07 3.00 

0.23E-09 3.00 

0.28E-05 ( 

0.15E-07 3.81 

0.23E-09 3.00 

0.36E-1 1 2.99 

CFL 

8 

4 

1 


Table IV 

L -2 Error of C-L method for different choices of r = 2u>oot with k = \,CFL = 1 


N 

a = 8 

a = 2 

a = 1 

a = 0.9 

a = 0.5 

16 

0.31E-05 

0.28E-05 

0.65E-05 

0.83E-05 

0.32E-03 

32 

0.15E-07 

0.15E-07 

0.15E-07 

0.15E-07 

0.18E-01 

64 

0.23E-09 

0.23E-09 

0.23E-09 

0.23E-06 

unstable 

128 

0.36E-11 

0.36E-1 1 

0.36E-1 1 

unstable 

unstable 
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From Table IV, we can see that for r < 2u; 0) 0-L becomes unstable while for r > 2u> 0 , 
the convergent of the scheme confirms the theoretical prediction. 


Case 2: Nonlinear scalar PDE 


Consider the scalar nonlinear initial boundary value hyperbolic equation 

Ut = U* —2irk cos(27 rk(x + /))(1 + sin(2ir&(x + /))) 

-1 < x < 1 , t. > 0 


(5.6) 


with initial condition 


U{x, 0) = 2 + sin(27rfcx) 
and boundary condition at x = 1 

f/( 1 , /) = g(t ) = 2 + sin(27rA:(l + t)) . 

This PDE has an exact solution given as U(x,t) = 2 + sin(27rfc(x + /)). 

Table V 

L ‘2 Error of C-L method for different choices of r = 2u;ocv with k — 1, cfl = 1 


N 

o = 8 

o = 4 

o = 3 

o = 2.5 

Exact BC 

16 

0.86E-02 

0.10E-01 

0.18E-01 

0.27E-01 

0.72E-02 

32 

0.40E-07 

0.40E-07 

0.40E-07 

0.11E-04 

0.39E-07 

64 

0.68E-09 

0.68 E-09 

0.68E-09 

unstable 

0.67E-09 

128 

0.11E-10 

0.11E-10 

0.11E-10 

unstable 

0.11E-10 


Different values of k are also tested, similar results are obtained. 
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